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Abstract 

In this paper, we present a method to measure the frequency and the frequency change rate of a 
digital signal. This method consists of three consecutive algorithms: frequency interpolation, phase 
differencing, and a third algorithm specifically designed and tested by the authors. The succession 
of these three algorithms allowed a 5 parts in 10^'^ resolution in frequency determination. The 
algorithm developed by the authors can be applied to a sampled scalar signal such that a model 
linking the harmonics of its main frequency to the underlying physical phenomenon is available. 
This method was developed in the framework of the Gravity Probe B (GP-B) mission. It was 
applied to the High Frequency (HF) component of GP-B's Superconducting QUantum Interference 
Device (SQUID) signal, whose main frequency fz is close to the spin frequency of the gyroscopes 
used in the experiment. A 30 nHz resolution in signal frequency and a 0.1 pHz/sec resolution in 
its decay rate were achieved out of a succession of 1.86 second-long stretches of signal sampled at 
2200 Hz. This paper describes the underlying theory of the frequency measurement method as well 
as its application to GP-B 's HF science signal. 
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I. BACKGROUND AND AVAILABLE SIGNAL 



The GP-B experiment aims at testing in Earth orbit two predictions of Einstein's general 
relativity using precision gyroscopes. This idea was independently proposed by Pugh |l| 
and Schiff [2I jsl in 1960. Both of these authors pointed out that according to the general 
theory of relativity, the angular momentum axis of a gyroscope in orbit about the Earth will 
precess about a direction normal to the orbital plane due to the gravitational interaction of 
the spinning gyroscope with its orbital motion, and simultaneously about the direction of 
the Earth's rotation axis due to the interaction of the spinning gyroscope with the angular 
momentum of the Earth. The first effect is known as the geodetic effect, and the second 
is known as the frame- dragging effect. On a 640 km polar orbit, the gyroscope drift rate 
due to the orbital motion about the Earth is 6.6 arcsec/yr (32 yurad/yr), while the orbital 
average drift rate due to the Earth's angular momentum is 0.041 arcsec/yr (0.20 /irad/yr). 

GP-B uses four gyroscopes spinning in a quasi torque-free environment and placed inside 
a drag-free satellite. The orientation of the gyroscopes is known thanks to an on-board 
telescope pointing towards a distant guide star, whose orientation with respect to an extra- 
galactic source is known. Furthermore, the satellite rolls about the telescope axis. As the 
orientation of the satellite with respect to an inertial reference frame is nonetheless known, 
reference frames linked to the satellite are qualified as inertial throughout this paper. GP- 
B's scientific goal can be fulfilled by measuring the orientation of the gyroscope's spin axis 
with respect to the satelhte. 

The gyroscopes are superconductive, which allows tracking of the orientation of their 
angular velocity vectors. Indeed, a spinning, superconducting body creates a magnetic 
dipole parallel to its spin axis, the London moment j4| js]. This magnetic dipole is then an 
excellent indicator of the direction of the instantaneous spin axis. Since the gyroscopes are 
almost perfectly spherical and uniform (A/ < 10~^), the spin axis direction is a very good 
indicator of the direction of the angular momentum. 

Low-noise Superconducting QUantum Interference Device (SQUID) magnetometers are 
thus used to measure the magnetic flux through a pick-up loop placed around each gyroscope 
created by the London moment plus a contribution due to a magnetic field trapped in the 
rotor. The SQUID signal is proportional to the magnetic flux through the pick-up loop. The 
SQUID output is an analog signal which, on board the satellite, is split into a low frequency 
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(LF) channel — which contains the London moment contribution — and a high frequency (HF) 
channel. Both HF and LF channels pass through a 780 Hz low pass analog filter. The LF 
channels then passes through an additional 4 Hz analog low pass filter and an additional 
gain stage. In this paper, we are only concerned about the HF channel, which is sampled 
at 2200 Hz and digitized with a 16 bit ADC with a range of ± 10 V. This ADC provides 
a resolution higher than the signal-to-noise ratio which is ~ 10^ near the gyroscope spin 
frequency. We refer to this digitized high resolution signal as the "HF SQUID signal". 

When each gyroscope transitioned below its critical temperature, the flux due to the 
residual magnetic field surrounding it before the transition was trapped on its surface, form- 
ing a large number of small magnetic sources. These sources are called fiuxons, and can be 
pictured as rigidly linked to the surface of the body, as shown for instance in |9|. A conse- 
quence is that the fiuxons exactly follow the motion of the body and create modulations, on 
the order of a few volts, at a frequency close to the spin frequency of the gyroscope: these 
modulations constitute the HF signal. 

The HF SQUID signal is sent intermittently in the form of 1.86 second long stretches 
that we refer to as "snapshots". Each snapshot contains 4096 points. We also have access 
to this signal in the form of Fast Fourier Transforms (FFTs) performed at regular intervals 
on the SQUID signal by the on-board CPU. The FFT is applied every 10 seconds to sets of 
4096 points of raw data, and a compacted form of its output is sent to the ground. 

An on-board FFT algorithm is applied to 1.86 second-long stretches of HF SQUID signal 
and can thus resolve frequencies 1/1.86 sec =0.54 Hz apart. This frequency resolution, 
called a 'bin', is poor: the FFT data is thus further processed. This additional processing is 
performed on the ground, and requires computing the FFT at a few frequencies. The value 
of the FFT at the first five harmonics of the signal frequency and the one corresponding 
to the 110 Hz calibration signal are thus telemetered. In addition, the value of the FFT at 
the two bins adjacent to each one of these six frequencies are also sent down, as well as the 
zero- frequency term. For each 1.86 second- long data stretch, the on-board FFT algorithm 
thus provides 19 data points. Since the high frequency data is sampled at 2200 Hz, the 
initial data stretch contains 4096 points and so does its FFT, so transmitting 19 data points 
is a dramatic reduction in bandwidth. 

Additionally, raw snapshots of the SQUID signal are also sent down on an irregular but 
frequent basis. These snapshots also consist of 1.86 second-long data stretches sampled at 
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FIG. 1: (Color online) Snapshot (time-series and spectrum) for Gyroscope 1, taken November 10, 
2004 



2200 Hz. An example of the time history of a portion a snapshot is provided in the upper 
pane Figure [H while the FFT of the entire shapshot is shown in the lower pane. 

In this paper, we describe the method that was used to determine the frequency of this 
HF signal, and therefore the rotor spin speed, with a resolution better than 30 nHz. This 
was critical for the GP-B data analysis because it allowed the time-varying readout scale 
factor to be determined to 1 part in 10^, the orientation of the spin axis with respect to the 
spacecraft to be determined to ~ 3 marcsec in 1 orbit, and the gyroscopes' relativistic drift 
rates to be determined to 20 marcsec/yr 6|. The time- varying scale factor is caused by the 
trapped flux contribution to the magnetic flux through the pick-up loop, which varies at the 
rotor spin ± spacecraft roll frequency, the rotor polhode frequency and at low frequency. 
With an accurate estimate of the rotor polhode and the estimate of the rotor spin speed to 
30 nHz (discussed here) the body-fixed orientation of the gyroscope rotor with respect to 
the spacecraft was determined with an accuracy of ~ 1 deg throughout the entire science 
mission, lasting 1 year. This information was necessary to determine the distribution of 
trapped flux on the surface of the rotor and the time varying readout scale factor [3]. 

Before delving into the 3 successive frequency estimation algorithms, we introduce no- 
tations specific to the GP-B experiment in order to explain the relationship between the 
gyroscope motion and the SQUID signal. We then describe how frequency interpolation 
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and phase differencing, a time-domain technique, were apphed to the FFT data to achieve a 
5 /iHz frequency resolution. We then show that this resuh was checked using the snapshot 
data. Variations of these first two techniques are known and their accuracy in determining 
monotone signals in the presence of Gaussian white noise and simple systematic effects have 
been studied iSl] . Finally, we show how the snapshot data and the result of the phase differ- 
encing were used to run the algorithm developed by the authors which allowed a 5 parts in 
10^" accuracy in frequency determination. 



II. GYROSCOPE MOTION AND HF SQUID SIGNAL FREQUENCY 

As the fluxons attached to the gyroscope move with it, HF modulations are created in 
the SQUID signal. The fluxons create a body-fixed distribution of potential: a model for 
this distribution can thus be written in body-fixed frame, for instance using a spherical 



harmonics expansion [10|, The SQUID measurement however takes place in the satellite 
frame - which, as explained above, is considered inertial. Therefore, if the appropriate set 
of Euler rotations is applied to rotate the body-fixed frame into the inertial frame, it is 
possible to express the HF SQUID signal as a function of the coefficients of the model of 
the magnetic potential distribution. 

The vectors Ji, I2 and are the principal inertia axes of the gyroscope and define an 
orthonormal body-fixed reference frame. We call (xj, yi, zi) an orthonormal inertial reference 
frame such that zi is aligned with the angular momentum L of the gyroscope. The angular 
momentum is inertially fixed as the gyroscopes are in torque-free motion. These notations 
are shown in figure [2l 

The first Euler rotation from (Ji, I2, 13) to (xj, yi, zi) is an azimuthal rotation by an angle 
(j)p. This rotation transforms (Ji, I2, 13) into a reference frame x'y'z' whose z' axis is aligned 
with I3. A polar rotation by an angle 7 is then applied: the new reference frame x"y"z" is 
such that its third axis z" is aligned with the angular momentum. The third Euler angle 
0s measures the angle by which the gyroscope has spun about L since a fixed time origin: 
a rotation of angle —0s about L is thus also needed to obtain the inertial reference frame 

(Xi,^!,^!). 

Consequently, the frequency at which any fluxon passes through the pick up loop is 
0s/27r -|- /p, where /p is the polhode frequency, such that /p ~ 0p/27r. As the modulations 
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FIG. 2: (Color online) Rotations from body fixed frame (Ii, 12,13) to inertial frame 

in the HF SQUID signal are due to the motion of the fluxons with respect to the pick-up 
loop, the frequency fz of this signal is thus given by: 



2n 



+ fp 



(1) 



An expression for (pg can be obtained from 



13 , formulas 89.4 and 90.4: 



L 



Ji 1 + a^sn2(r, k'^) ' 
where the characteristic, a, and the elliptic modulus, k, only depend on the moments of 
inertia Ji, I2 and Is, the rescaled time,r, depends on the three moments of inertia and the 
angular momentum, L (see 10|) and sn is the Jacobi elliptic integral referred to as "sinus 
amplitude" and defined by: 



1 



sn(r, k) = X when t = [ — 

Jo J(l-s^)(l 



^0 V(l -s2)(l- A;V) 

The first term in Eq. ([2]) is the larger by a factor of 10^, and the second term is a modulation 
at twice the polhode frequency. The three Euler rotations are summarized in the definition 
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of the spin axis uj: 

^ = <i>ph + iy' - 4>sz". 



The spin frequency fs is thus 



(27r/,)2 = 02 + 02 ^ ^2 ^ 20,0pcos7. (3) 



As (j)p and 7 can be as large as fp which is on the order of 0.1 niHz on all gyroscopes [10|, 
(ps lies within 0.1 mHz of the spin frequency fs and is thus on the order of 100 Hz. The HF 
SQUID signal frequency fz is thus close to the spin frequency /«. 

The first two steps of the frequency determination method presented in this paper are 
applied to the HF SQUID signal in order to determine fz- The third step of the frequency 
determination method presented in this paper thus aims at determining 0^ with a 20 nHz 
resolution. In [l^, a procedure is given to measure the polhode frequency fp with a 10 nHz 
accuracy. Note that this represents an accuracy of 1 part in 10^ as the polhode frequency 
is on the order of 0.1 mHz. Therefore, from ([1]), an absolute accuracy better than 30 nHz 
in the the determination of the frequency of the HF SQUID signal can be achieved. In 
relative terms, the claimed accuracy of the frequency fz determination is thus on the order 
of 5 part in 10^°. 



III. MILLIHERTZ LEVEL SIGNAL FREQUENCY DETERMINATION 
THROUGH FREQUENCY INTERPOLATION 

By interpolating the FFT data, it is possible to go beyond the 0.54 Hz frequency resolution 
and to determine the frequency with a 1 mHz precision. This method is well known and 
has been used for several decades {ll]. It consists of approximating the signal by a pure 
sine around its main frequency, and finding the frequency of such a sine if it had the same 
amplitude diagram as the FFT of the signal. We now give a mathematical description of 
this method. 

As the interpolation is applied to the FFT data which is itself computed using short 
stretches of SQUID signal, we neglect the polhode harmonics in this description. We show, 
without loss of generality, how frequency interpolation is implemented for a single frequency 
signal. 
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Let ^hf(^) be our signal, so that 



zhfW = As cos (27r/^t + 



Its discrete Fourier transform is: 

Aft -r J, 



nf) 



-e 



As 



^ ■^_^-iS<t>^-i27r(f+f,)At^ 



sin[27r(/ - f,)AtN/2] 
sin[27r(/ - f,)At/2] 
sin [27r(/ + /,)AtiV/2] 
sin [27r(/ + /,)At/2] 



(4) 



where: At is the time between two consecutive samples (1/2200 sec for a 2200 Hz sampling 
rate), N = 4096 is the number of points in the FFT, NAt is then the total snapshot duration, 
fd = l/{NAt) is the bin frequency, / is the frequency at which the FFT is computed. As is 
the signal amplitude and 6(j) is the phase shift. 

Near the HF signal's frequency, / ~ fz, so the first term in (jl]) is dominant and we neglect 
the second one. Since the FFT is discrete, we only obtain its values at multiples of the bin 
frequency fa- Furthermore, as mentioned above, we have at our disposal the value of the 
FFT at the central bin and at the two adjacent bins. Let's define the integer n such that 
nfd is the multiple of fd closest to the signal frequency f^. We note: 



Fn = F{nfd) Fn+i = F[{n + l)fd] = - 1)/^]. 

Then the quantity we use in the interpolation is: 



Sin [27i{nfd - /.) AtiV/2] 
sin [2n{nfd - /.) At/2] 

and F^+i are defined similarly. Let's also introduce the ratios: 



7„ = ^e'^<^e 
2 



Rn+l 



n+1 



Rn-1 



The values of those two ratios are obtained from measurements. We define the quantity Xn 
to be, 

= 27f{nfd - h)NAt/2. (5) 

Then, by performing a Taylor series expansion of Fm in the quantity Xm/N, where m takes 
values n — 1, n, + 1, we obtain. 



Rn 



+1 



Xn + TC 



0(f) fl..- 



Xr,. - vr 



(6) 



As implied by Eq. [5l the term (xn/N) here is 10~^ at most, since near the main peak the 
signal frequency, fz, is at worst one frequency bin away from nfd- From these relations, we 
derive two formulas for the frequency fz of the HF SQUID signal, depending on the sign of 
Xn- For Xn > 0: 



Averaging the two formulas for fz we and obtain: 



fz = nfd + 



2NAt 

For Xn < 0: 



Rn+l Rn-1 



Rn+l — 1 -Rn-1 + 1 



+ O (10-3) (9) 



Averaging these two formulas gives: 

+ C(lO-3) (12) 



fz = nfd + 



2NAt 



Rn+l Rn-1 



1 + Rn+l Rn~l — 

Therefore, the frequency interpolation yields a determination of the HF signal's frequency 
fz, whose error is 0{xn/N) < vr/A^, or 10^^ Hz. Indeed, the frequency computed by this 
method typically showed a 100 /iHz spread, as shown on figure [31 

The result derived in ( !T2|) was obtained for the simple case of a signal with only one 
harmonic. However, the HF signal contains many harmonics of its main frequency, and 
as explained in [l^, [l^, the odd harmonics have a significantly larger amplitude. The 
interpolation procedure is therefore separately applied to the fundamental frequency, fz, as 
well as to the 3'"'^ and harmonics of the FFT data. The results from each of those 3 
interpolations are then averaged in order to obtain a determination of the signal frequency 
fz- Results for a 6- hour long data stretch for gyroscope 1 are shown in figure [31 

The formula ([T2|) for the signal's frequency is valid when no window is applied. Applying 
windows to a signal consists of multiplying it by a carefully chosen function -for example, 
a sine square- so the frequency components away from the main frequency are attenuated 
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FIG. 3: (Color online) Frequency fz from interpolation, gyroscope 1, Aug. 19 2004 
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16| . In the case of the windowed signal, the analytical expression of fz is more complex 



but the principle of the derivation remains exactly the same. 



IV. MICROHERTZ LEVEL SIGNAL FREQUENCY DETERMINATION 
THROUGH PHASE DIFFERENCING 

The interpolation method relies on the use of the amplitude of each independent FFT 
and, in practice, yields a 1 mHz or better frequency resolution. We now present an additional 
method called phase differencing that uses the phase of the peak FFT bin, which is the one 
closest to the frequency of interest, at two different times. This method yields more than two 
additional orders of magnitude in frequency resolution. Phase differencing consists of using 
a precise measurement of the change in the phase of the FFT computed on two distinct 1.86 
second-long data stretches. Knowing the time elapsed between the two data stretches, we 
can determine the signal frequency. However, to resolve the 2n ambiguity, an initial estimate 
of this frequency is necessary. Indeed, let 0i be the phase of the FFT obtained at time ti 
and 02 the phase of the FFT obtained at time ^2- We can assume that the frequency fz of 
the signal is constant if ti and ^2 are close (we indeed show later that the characteristic time 
of the decay in fz is 7,000 - 25,700 years, depending on the gyroscope). We then have: 



11 



02-0i = 27r/,(t2-ti)-27rfr, 



(13) 



where K is an integer such that 02 — 4>i is between and 27r. We then need to know the 
value of K in order to extract the frequency from the previous equation. This requires 
knowing an estimate /|** of the signal frequency so that: 



The value of K then has to be chosen to ensure condition f|T4l) . In order for K to be known, 
/l'^* must lie within 1/10 Hz of the real value / since {t2 — ti) = 10 sec. This resolution is 
easily achieved by the interpolation algorithm. Furthermore, this 1/10 Hz requirement is 
more stringent than the 0.54 Hz precision of the FFT: the result of the interpolation is thus 
needed to carry out the phase differencing. 

Then, once K is computed, equation f|T3|) can trivially be solved for the signal frequency 
fz- This new estimate has a very good precision. Indeed, the phase is known with a 
conservative error of 10 arcsec (5 x 10^^ rad), which at the gyroscope spin speed corresponds 
to a timing error of approximately 7 x 10~^ sec, which may be compared to the short-term 
on-board clock accuracy of approximately 10~^ sec. The main error term is thus the phase 
measurement, and for a {t2 — ti) = 10 sec interval, we obtain a precision better than 5 /xHz. 
Figure H] shows the extra resolution that can be achieved with phase differencing. It should 
be compared with figure [3l which only shows the output of the interpolation on the same 
data stretch. 

The large modulations in fz are at polhode frequency and are a systematic error in the 
estimation method. Indeed, we used a simple model for the HF signal zhf which only takes 
into account the harmonics of fz, whereas harmonics of the polhode frequency also contribute 
to the signal's spectrum. These frequencies then leak into our current measurement of fz- 
The model introduced in the last section of this article takes into account all the frequency 
components of the HF signal, and yields an estimate of fz free of large variations at polhode 
frequency. 

As a conclusion, by interpolating each FFT and using phase differencing between FFTs, 
we can obtain an estimate of the signal frequency with an error of a few /xHz. We have thus 
gained 5 orders of magnitude as compared to the FFT bin width of 0.54 Hz. Furthermore, 
this determination only relies on estimates of the FFT for a small number of frequencies 
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(14) 
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FIG. 4: (Color online) Frequency fz from phase differencing, gyroscope 1, Aug. 19 2004 



gyro signal frequency fz (Hz) 

1 79.38715 

2 61.81759 

3 82.09202 

4 64.85030 



TABLE I: Typical values for signal frequency fz obtained from phase differencing 

being transmitted from the spacecraft. It then has the double advantage of requiring little 
bandwidth and using less on-board CPU time, since the interpolation and the phase differ- 
encing are carried out on the ground. Typical values for fz obtained after phase differencing 
are given in table [H They have been measured on Feb. 6*'* 2005 at 07:39 GMT. 

Figure [5] shows a linear decay of the frequency fz for gyroscope 4 over four months. 
This decay was observed on all four gyroscopes. As the polhode frequency fp increased 
slowly throughout the experiment, equation ([T]) implies that the decay in frequency can be 
attributed to a decay in 0s- 

The FFT is taken on-board every 10 seconds for approximately 40% of each orbit, a 
period called Guide Star Invalid (GSI). During the rest of the orbit, the Guide Star Valid 
(GSV) period, the on-board FFT algorithm is turned off in order to allow the CPU capacity 
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FIG. 5: (Color online) HF signal frequency fz on gyroscope 4 over 4 months 

to be used to track the guide star: we then have about 4000 frequency measurements per 
day with a 10-second spacing, each one with a 5 yuHz standard deviation. A covariance 
analysis yields a 0.1 nHz/sec uncertainty on the value of the signal frequency decay rate 
using 24 hours worth of data (assuming 16 GSI periods per day, each lasting 40 minutes and 
containing 240 samples, and that for each GSI period a frequency offset must be estimated 
together with the decay rate). An averaging time of 24 hours is chosen because there are 
sometimes large gaps in the FFT data, ranging from a few hours to a few days, that occurred 
at roughly 1-2 day intervals that preclude longer averaging times. These gaps are related to 
details of the spacecraft operations and, on rare occasions, to spacecraft anomalies. 

Typical values for the frequency decay rate and the characteristic time of the decay are 
given m table HI 

V. VERIFICATION USING SNAPSHOT DATA 

The results obtained so far are derived from the analysis of the FFT data. To increase the 
level of confidence in these results, they have been checked by an independent verification 
using the snapshot data. These data are more voluminous, hence much richer, than the 
FFT data: a snapshot comprises 4096 points, whereas an FFT compacts it into 19 points. 
Typically, snapshots are transmitted about every 40 seconds during approximately an hour. 
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gyro 


frequency decay rate 
(nHz/sec) 


characteristic time 
of decay (years) 


1 


0.16 


15,800 


2 


0.14 


13,400 


3 


0.36 


7,000 


4 


0.08 


25,700 



TABLE II: Typical decay rate and characteristic time for the frequency fz of the HF signal 

separated by gaps ranging from an hour to up to two days. 

A two-step algorithm similar to the procedure followed in the FFT analysis was imple- 
mented in order to determine the signal frequency from the snapshot data. First, a Fast 
Fourier Transform of each snapshot was performed and interpolated in order to produce an 
initial estimate of the spin frequency. A two-point interpolation algorithm was used (instead 
of three-point) to simplify the coding. A least-squares fit of each snapshot to a sum of up 
to 51 harmonics of the frequency obtained from the interpolation was then performed. 

Secondly, a nonlinear fit to the snapshot data for the signal frequency, fz, and the am- 
plitudes of the sine and cosine components of 27cfz t and its harmonics (up to the 51st) 
was performed using the NonLinearRegress routine available in Mathematica. The initial 
conditions were, for fz, the value obtained after the interpolation and, for the harmonics 
coefficients, the output of the linear fit. 

This two-step analysis yielded an independent estimation of the signal frequency as well 
as a large number of harmonics coefficients using a data set and a method which was different 
from those used in the FFT analysis. 

Figure [6] compares fz determined from the FFT and snapshot data on gyroscope 4 for 
February 7*^ 2005. Note that, like for the FFT data, the estimation of fz obtained from 
the snapshot data shows variations at the polhode frequency. This is expected as the model 
for the signal used in this determination also neglects the harmonics of polhode. These 
results are typical and have been observed on all the gyroscopes throughout the mission. 
The signal's frequency obtained from the snapshots is typically noisier than the one obtained 
from the FFT (the spread is on the order of 30 /xHz). The higher level of noise is likely due 
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FIG. 6: (Color online) fz obtained from FFT and snapshots, gyroscope 4, Feb. 7 2005 

to the on-board analog-to-digital converter, which exhibited noise at all multiples of 10 Hz. 
Nevertheless, both time histories coincide strikingly despite the difference in data sets and 
in algorithms, which increases our confidence in the measurement. 

Since more information is contained in the snapshot data, including harmonics higher 
than the 5**^, it was used as a reference for the remainder of the analyses presented in this 
paper. The analog-to-digital converter noise is later estimated and subtracted from the 
snapshot data to remove its effects. A more detailed discussion of the justification for this 
choice can be found in 10|. The last step of the method indeed dramatically increases 
the signal's frequency resolution, which fully offsets the poorer quality of the determination 
obtained from the snapshots. 

VI. NANOHERTZ LEVEL SIGNAL FREQUENCY DETERMINATION 
A. Overview 

The interpolation and phase differencing methods allowed a determination of the signal 
frequency fz with a microhertz level accuracy. As the error in the determination of fp is 10 
nHz, formula ([1]) indicates that a better resolution in fz can be achieved if (pg is measured 
accurately. 

In section HV] of this paper, it is shown that the linear decay in the measured frequency 
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fz is due to a decay in (pg. The expression given in ([2]) for (pg shows that the decaying term 
is necessarily the angular momentum L. We can then write: 

/3-/1 1 



0, = 27r(Ci + C2(t-to)) 



1 + 



(15) 



where to is a chosen time origin. 

An accurate determination of the coefficients Ci and C2 is thus necessary to estimate 0^. 
This last section describes how this estimation was carried out. The underlying principle is 
to expand the HF SQUID signal in harmonics of 0=. A model for the complex amplitude 



12] and is given in its compact form in 



Hn{t) of the n*^ harmonic of (ps has been derived in 
equation (|T71) : the values Ci and C2 yielding the time history for 0s (t) such that the measured 
harmonics Hn{t) are best accounted for by the model are taken as the best estimates of Ci 
and C2. 

Once Ci and C2 are known, the value of (ps can be computed from (fT5l) . and the signal's 
frequency fz can be obtained from ([1]). 



B. A cost function to estimate 



Reference gives an expansion of the HF SQUID signal zhf in harmonics of 0^: 

00 

^HF(t)= ^nWe-^^^-^^W), (16) 

n=— oo,n^0 

A model for Hn(t) with n odd can be derived using a similar procedure to that used in 
[l^ . and this has been done in 1^, Q. This model states that Hn{t) depends on time only 
through the angles 0p(t) and 7(t), determined in jlij], so that: 

i/„(t) = i/„(0p(t),7(t)). 

This model is furthermore linear in the coefficients of the expansion of the magnetic po- 
tential in the body-fixed frame. In other words, -?/„(0p(t), 7(t)) is linear in a state vector 
A whose coefficients are constant as they only depend on the body-fixed magnetic poten- 
tial distribution. We call M„(0p(t), 7(t)) the observability matrix for this linear model. Its 
expression is given in |l^. We can thus write: 

i^n(0p(t),7(i)) = M„(0,(t),7(t))A (17) 
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Therefore, if the Euler angle (psit) is known accurately, Hn{t) is known from (fT6ll . and by 
running a linear least squares algorithm, we can estimate the vector A that minimizes the 
residuals Ji where: 



Ji 



mm 

A 



(18) 



-M„(7(t),0p(t))A 

However, errors in the estimation of Ci and C2 lead to an approximate value (j)1^^{t) for 
the third Euler angle. The coefficient H^'P{t) obtained from the expansion f|T6|) of the HF 
signal in harmonics of (f)g'^^{t) is therefore also approximate. Consequently, the model (fT7|) 
does not apply exactly: the vector A"'^^ that solves the least squares problem ( |T71) with the 
approximate value H'^'P{t) for Hn{t) yields a minimum norm Jf^^ of the residuals which is 
larger than the value Ji defined in ( fT8|) . 

Therefore, the best estimate 0^ of the third Euler angle is the one that yields the smallest 
value for the residual J^^^ . In other words, 0^ minimizes the functional J defined by 



mm 

A 



/f„[0,](t) -M„(t)l 



(19) 



The notation J 



means that the argument of the functional J is the function 



In reference lOj, an expression for 0^ is given based on the integration of the expression 
ISD for 0,: 



At) 



s, + 27rCi [{t - to) + ni(t, to, /i, I2, /3, fp)] 



-271C2 



-{t-to)^ + U2{t,to,h,l2,hJp) 



(20) 



In this expression, 0^^ is 0^ at the time origin to- The terms Hi and 112 are integrals of 
the Jacobi elliptic function sn and depend only on the moments of inertia and the polhode 
frequency. The function 0^ is thus fully determined by three variables: 0^0, Ci and C2. 

Using equation (120|) . we see however that the initial value 0sq is just a constant shift to 
the angle 0s (t). From the expansion f|T6|) . a constant shift in 0s (t) yields a constant phase 
shift e*"'^*o in Hn[(f)s]- In the expression (fTOll . we can factor out this constant phase shift, so 
that the term A becomes Ae'^^^^o . And since J[0s] is a minimum over all values of A, its 
value is then independent of (pso- Consequently, the best estimates of Ci and C2 minimize 
J[0,]. 

In other words, the estimation of Ci and C2 is a minimization problem of a function J2 
of two variables: 



J2(Ci, C2) = min i7„(Ci, C2, t) - M^t), M^))^ 

A 



(21) 
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We have thus constructed a cost function J2 which can be minimized in order to estimate 
the two coefficients Ci and C2. 

This frequency determination method is general. Indeed, it only supposes that a model 
for the harmonics of the signal's frequency is available: the frequency yielding the harmonics 
best ffi by the model is then the best estimate of the signal's frequency. 

Note that in the specific case of the GP-B gyroscopes, the modeled quantities are not 
exactly the harmonics of the signal's frequency fz of the signal but rather of (ps = '^'^{fz — fp), 
where fp is known accurately. 

C. Implementation 

We now introduce the algorithm developed by the authors to estimate Ci and C2. 

A nonlinear minimization routine based on a modification of the Nelder-Mead (NM) 

n 

simplex p/7|, used in Matlab's fnimsearch.m function, was implemented. A simplex is a 
set of n + 1 points in an n-dimensional space such that the points are not degenerate. In 
2-space the^ simplex is a triangle. A detailed account of how this algorithm works can be 
found in [14]. This nonlinear simplex routine was chosen because of its high accuracy for non- 
smooth objective functions compared with gradient or Gauss-Newton methods, implemented 
in Matlab's Isqnonlin.m function for example. The Nelder Mead simplex is a nonlinear 
minimization algorithm and as such requires the input of initial conditions. Since we know 
a priori the region of possible values for Ci and C2, we utilize a NM simplex method that 
is modified to include a constraint on the search domain. From equation f fTSj) . we find that, 
for the GP-B gyroscopes: 

0,(to) = 27rCi(l + 0(10-6)) 

And since: 

05 = 27r(/^ - fp) 

we find: 

Ci = (/.(to)-/p(to))(l + 0(lO-'^)) 

Using the measurement of fz obtained previously, we thus have an initial condition for Ci 
as well as a search interval whose width is on the order of 100 /iHz. 
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Similarly, we use the measurement of the decay in to find an initial condition and a 
search interval for C2. As \dfp/dt\ can be as large as 1/10 dfz/dt, the relative size of the 
search interval is larger than for the coefficient Ci. we opted for a search interval whose size 
is 25% of the measured decay rate of f^. 

Using these two initial conditions and search intervals, we define a search region in the 
(C'ijC'2) plane. The NM simplex algorithm will then search for the values of Ci and C2 
within this region that minimize the cost function J2. 

The NM algorithm is implemented in Matlab's built-in fminsearch.m routine. The al- 
gorithm first evaluates the cost function J2 on a set of three points A, B and C, called a 
simplex, within the search region. The Matlab algorithm finds on which vertex A of the 
simplex the function has the maximum value. Let B and C be the other two vertices of the 
simplex and H the center of segment BC . The algorithm evaluates the function at a point 
A' such that A A' is orthogonal to BC and A A' = 6 AH with 6 = 3. (see figure [ZD. This 
procedure is called an 'expansion'. Indeed, the simplex A' BC has a bigger area than ABC. 
If the cost function J2 has a lower value in A' than in A, the algorithm is searching in the 
right direction in the (Ci,C2) plane, and the procedure repeats, now starting from A'. If 
the function has a higher value in A' than in A, the expansion did not decrease the value of 
the function and should then not be pursued. The cost function is then evaluated at a new 
point A" inside the ABC simplex (see figure [8]), such that AH = pHA", where p = 2, and 
the procedure repeats. Matlab's algorithm stops when the dimension of the simplex, that is 
the length of its longest side, is lower than the specified tolerance. 

There are two main limitations in the available Matlab routine. First, only one tolerance 
can be specified. This is in our case a real issue as both the absolute and relative tolerance 
on Ci and C2 differ by many orders of magnitude. Furthermore, the algorithm does not 
provide any way to decrease the likelihood of honing in on a local minimum. 

Matlab's fminsearch.m routine was thus modified to address these two issues. First, a 
simple but crucial change was brought to fminsearch.m so that the termination tolerances 
on the two dimensions were allowed to differ: the modified NM algorithm terminates only 
when the length of the simplex along both the Ci and C2 dimensions are smaller than 
values yielding the required precision in (pg-. these numbers are referred to as "termination 
tolerances" . 

After examining several algorithms for decreasing the likelihood of reaching a local min- 
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FIG. 7: (Color online) Expansion of the simplex 

imum, the following steps were determined to be the most effective for finding the global 
minimum, when applying the above model to the snapshot data. The first one is to add an 
extra expansion when the original algorithm terminates: in case a value J2 lower than the 
current minimum is found, the algorithm restarts. Once a minimum is attained, a 20 x 20 
grid in the (Ci,C2) plane is then created, each point being spaced by the termination tol- 
erance. J2 is evaluated at each of these points. If one of the 400 points in this local grid 
yields a lower value of J2 than the current minimum, the algorithm is restarted at this point. 
The algorithm then terminates and returns the values of Ci and C2 that minimize J2. As a 
third and last step, we then restart this algorithm twice such that at each iteration we use 
the same initial search region and the current best guess of Ci and C2 as initial conditions. 
The new best estimates of Ci and C2 are then compared with the previous ones and it has 
been observed that 3 iterations are sufficient to obtain estimates within a few times the 
termination tolerances. 

The estimation of Ci and C2 was performed using the modified NM algorithm to minimize 
the cost function J2. This cost function is however based on a spherical harmonics model for 
the fiuxon distribution around the gyroscope, and the order of this expansion, conventionally 
denoted by the letter /, is in theory a free parameter. Nevertheless, early estimations and 
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FIG. 8: (Color online) Simplex honing on a minimum : case ^2(^0 > J2{A) 

fits to the SQUID data suggested a range of values for /. The NM search was thus repeated 
for all reasonable values of /. The optimal values for Ci and C2 were obtained by averaging 
the estimates corresponding to different values of /. 

The standard deviation of the ensemble of estimates of Ci and C2 was on the order of 
10 nHz for Ci and 10~^ nHz/sec for C2. These results were obtained using one-day batches. 
Therefore, a conservative estimate of the uncertainty in 0^ at the end of a one-day batch is: 
10 nHz + 10~^ nHz/sec x 86,400 sec ~ 20 nHz (see Eq. (|T5l) ). The uncertainty in the signal 
frequency, = 0s/27r + fp, is then conservatively estimated to be 30 nHz, for a 10 nHz 
uncertainty in the polhode frequency, fp. This uncertainty bounds the uncertainties for all 
gyroscopes during the entire science mission. For a gyroscope spinning at 60-80 Hz, this 
represents a relative error of 5 x 10^^°. Unlike the previous two approaches, the polhode 
motion of the gyroscopes is modeled explicitly. Therefore, the polhode frequency harmonics 
are completely separated from the gyroscope spin frequency estimates. The value of the 
signal frequency and decay rate for February 6*^ 2005 are given in tables IIIII and IIVI They 
can be compared with the results obtained with the phase differencing method on the same 
day in tables [T] and [Tll 
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gyro signal frequency (Hz) 

1 79.38746144 

2 61.81765160 

3 82.09121932 

4 64.85036174 



TABLE III: Signal frequency on February 6*'^ 2005 at 07:39:00 GMT 

gyro frequency decay rate (nHz/sec) 

1 0.1581 

2 0.1432 

3 0.3634 

4 0.0803 



TABLE IV: Signal frequency fz decay rate on February 6*" 2005 
VII. CONCLUSION 

This paper describes three successive algorithms for estimating the frequency of a digital 
signal. The first two of them, interpolated FFT and phase differencing, known previously, 
were properly adjusted for the needs of the analysis of the HF signal obtained in the Gravity 
Probe B experiment. The third method is new; it is based on the accuracy achieved by the 
previous two approaches, and involves a nonlinear estimation of a slowly changing frequency 
when many harmonics of another time- varying frequency are present. The last method 
requires accurately modeling the rigid body motion of a nearly torque-free gyroscope, and 
utihzes the specially modified Nelder Mead simplex method. 

These algorithms are applied to a particular problem of precisely determining the spin 
frequency of the ultra-precise GP-B gyroscopes. The measured HF SQUID signal was used 
for this analysis whose signal to noise ratio was ~ 10^. We achieved the ultimate relative 
accuracy of 5 parts in 10^°, corresponding to the 30 nHz absolute accuracy for a gyro spinning 
at 60-80 Hz. This accuracy level was necessary for the estimation of trapped magnetic flux 
on the surface of each gyroscope rotor to 1%, the resulting time variations in the gyro readout 
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scale factor to 10 ^, and ultimately the relativistic gyro drift rate to 20 marcsec/yr. 
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